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ABSTRACT 

We present self-consistent numerical simulations of the sun's convection zone 
and radiative interior using a two-dimensional model of its equatorial plane. The 
background reference state is a one-dimensional solar structure model. Turbulent 
convection in the outer convection zone continually excites gravity waves which 
propagate throughout the stable radiative interior and deposit their angular mo- 
mentum. We find that angular velocity variations in the tachocline are driven by 
angular momentum transported by overshooting convective plumes rather than 
the nonlinear interaction of waves. The mean flow in the tachocline is time de- 
pendent but not oscillatory in direction. Since the forcing in this shallow region 
can not be described by simple linear waves, it is unlikely that the interaction of 
such waves is responsible for the solar cycle or the 1.3 year oscillation. However, 
in the deep radiative interior, the interaction of low amplitude gravity waves, 
continually excited by the overshooting plumes, is responsible for the angular 
velocity deviations observed there. Near the center of the model sun the angular 
velocity deviation is about two orders of magnitude greater than that in the bulk 
of the radiative region and reverses its direction (prograde to retrograde or vice 
versa) in the opposite sense of the angular velocity deviations that occur in the 
tachocline. Our simulations thus demonstrate how angular velocity variations in 
the solar core are linked to those in the tachocline, which themselves are driven 
by convective overshooting. 

Subject headings: convection: overshoot, mixing, internal gravity waves, solar 
interior 
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1. Introduction 

Internal gravity waves are ubiquitous in nature. Their influence can be observed in 
striated cloud structures in our own atmosphere many days of the year. In addition to 
the visual display of these waves in our atmosphere, they have several other more profound 
consequences. The interaction of gravity waves with angular velocity shear produces the 
quasi-biennial oscillation (QBO), which dominates variability in the equatorial stratosphere 
and affects ozone levels. 

Gravity waves have been invoked to explain several physical problems in the sun's in- 
terior such as: (1) providing the extra mixing required to solve the Li depletion problem 
(Garcia-Lopez & Spruit 1991) (2) increasing turbulent mixing and affecting the solar neu- 
trino production (Press 1981), (3) maintaining the solid body rotation of the sun's radiative 
interior (Schatzman 1993, Kumar & Quataert 1997) (4) and controlling the solar cycle (Ku- 
mar, Talon & Zahn 1999). They have also been used to explain the orbital properties of 
binary star systems (Terquem et al. 1998). These theories, however, are only as good as the 
poorly understood excitation and evolution of gravity waves in stellar interiors. What grav- 
ity wave spectra and amplitudes are generated? How effectively do gravity waves transport 
angular momentum and mix species? 

Some confusion has existed. For example, although gravity waves were postulated to 
extract angular momentum from the solar radiative interior and so enforce solid body rotation 
(Schatzman 1993, Kumar & Quataert 1997), it was quickly pointed out (Gough & Mclntyre 
1998, Ringot 1998) that gravity waves tend to enhance local shear rather than smooth it. 

Recognizing this ant i- diffusive nature of gravity waves in shear flows, several authors 
published papers postulating an angular velocity oscillation at the base of the solar convection 
zone analogous to the quasi-biennial oscillation (QBO) in the Earth's stratosphere (Baldwin 
et al. 2001). In one paper (Kumar, Talon & Zahn 1999), a spectrum of gravity waves 
generated by the overlying convection is prescribed (Goldreich et al. 1994). This spectrum 
is then integrated to give a flux of angular momentum that is transferred from the waves 
to the mean flow, resulting in a periodic oscillation of angular velocity at the base of the 
convection zone with a timescale of about 20 years. In another paper (Kim & MacGregor 
2001), a two wave model is assumed: one prograde propagating wave, with a prescribed 
angular momentum flux, and one retrograde wave with another (negative) prescribed flux. 
In addition, these authors include viscous dissipation in the evolution equation for the mean 
flow and find that the nature of the resulting angular velocity oscillation depends sensitively 
on the value of the assumed viscous diffusivity. For large viscous diffusivity, a steady solution 
is found; whereas for small values a chaotic oscillation is found. A periodic solution is 
recovered only for intermediate values of the viscous diffusivity. 
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More recently, the Kumar, Talon & Zahn (1999) model has been extended to show that 
the oscillating shear layer at the base of the convection zone acts as a filter on the low fre- 
quency, short wavelength waves (Talon, Kumar & Zahn 2002, Talon & Charbonnel 2005). 
In their theory, the filter preferentially damps prograde waves, allowing predominantly ret- 
rograde waves, which carry negative angular momentum, into the deep interior. When these 
waves dissipate they transfer their negative angular momentum to the flow and therefore 
effectively extract prograde angular momentum from the low-latitude deep solar interior, 
leading to solid body rotation in the solar radiative zone. 

The major shortcoming of both of these models is twofold. First, neither model self- 
consistently calculates the generation of the gravity waves by the solar convection zone. In 
Kumar, Talon & Zahn (1999) a spectrum and amplitude of gravity waves is assumed, a 
spectrum that is unfortunately, untestable with observations and has not been reproduced 
in numerical simulations. In addition, this spectrum neglects the main production of gravity 
waves due to overshooting plumes. Second, and perhaps more importantly, neither model 
calculates the nonlinear wave-wave interactions which provide this "flux" of angular mo- 
mentum. Rather than parameterizing the nonlinear interaction of waves as a prescribed 
flux (Kumar, Talon & Zahn 1999 and subsequent papers, Kim & MacGregor 2001), this 
interaction should be self-consistently calculated. 

Here we present numerical simulations that address both of these issues. Our model 
solves the fully nonlinear Navier-Stokes equations in both the convective and radiative re- 
gions. Therefore, the gravity waves are self-consistently generated by an overlying convection 
zone and the nonlinear terms are retained in the radiation zone to account for the nonlinear 
wave interactions which affect the mean flow. 



2. Wave - shear flow interactions 

In this section we briefly review the hydrodynamic process by which waves can transport 
their angular momentum to the flow in which they travel. The natural example of this process 
is the QBO. For a more complete discussion of the QBO see, for example, Baldwin et al. 
2001, Lindzen (1990) and Holton (1994). 

For simplicity consider the Boussinesq equations for a viscous fluid in two dimensional 
cartesian (x-z) coordinates, x being the horizontal coordinate and z the vertical coordinate. 
After decomposing the horizontal velocity into mean flow (U, a function of z) and fluctu- 
ating (u') components and taking the proper horizontal average, the equation for the mean 
horizontal velocity becomes: 
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In (1), v! represents the fluctuating horizontal velocity and w' represents the fluctuating 
vertical velocity. This equation represents how small scale wave-wave interactions (Reynolds 
stresses) influence the mean flow (U). The vertically propagating waves carry horizontal 
momentum vertically. Eliassen & Palm (1961) showed that (in this simplified model with 
no viscosity) the first term on the right hand side of (1) is zero unless there is some wave 
attenuation (such as radiative damping or critical layers). Therefore, in the absence of wave 
attenuation, there is no transfer of angular momentum between the mean flow, U, and the 
Reynolds stresses. In the sun the proposed mechanism for wave attenuation is via radiative 
diffusion, while in the Earth the process is likely wave breaking. 

In the presence of differential rotation the picture becomes more complicated. In a 
rotating fluid, waves generated at a particular frequency are doppler shifted away from that 
frequency according to the equation: 



where m is the horizontal wave number, uj gen is the frequency at which the wave is generated 
in the frame rotating at Q gen and u>(r) is the frequency measured relative to the local rotation 
rate, Q(r). 

Where Q(r) > VL gen prograde waves (m > 0) are shifted to lower frequencies and retro- 
grade waves (m < 0) are shifted to higher frequencies. Since radiative damping is strongly 
frequency dependent (damping length oc cj 4 ), prograde (lower frequency) waves are damped 
in a shorter distance than retrograde (higher frequency) waves. This differential damping 
causes prograde waves to deposit their (positive) angular momentum closer to the genera- 
tion site than retrograde waves deposit their (negative) angular momentum, thus leading to 
a shear layer. As prograde (retrograde) waves continually deposit their positive (negative) 
angular momentum, the angular velocity amplitude and gradient continually increase. This 
increased angular velocity causes prograde waves to be shifted to ever smaller frequencies 
which are damped even closer to the generation site. In this way the peak in the prograde 
layer moves toward the source of the waves. When the prograde shear becomes sufficiently 
steep, it is broken down by viscous diffusion, leaving behind the retrograde layer. This pro- 
cess repeats, with the period of the process inversely proportional to the wave forcing. The 
prominent features of this physical process are then twofold: (1) prograde flow lies above 
retrograde motion (or vice versa) and (2) this pattern propagates toward the generation site. 

This physical picture was proposed initially by Lindzen & Holton (1968) and Holton 
& Lindzen (1972) to explain the QBO observed in the Earth's atmosphere. Later, the 
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physical theory was tested in the remarkable experiment by Plumb & McEwan 1978 and 
the basic physical mechanism was recovered. It has also been suggested that an oscillation 
similar to the QBO occurs in Jupiter's atmosphere (Leovy et al. 1991) and is coined the 
quasi-quadrennial oscillation (QQO) because of its four year period. 

The robustness of this mechanism led astronomers to hypothesize that the same physical 
mechanism could be operating in the solar tachocline (Kumar, Talon & Zahn 1999, Kim & 
MacGregor 2001, Talon & Charbonnel 2005). This could then provide a handsome explana- 
tion for oscillations at the base of the convection zone, whether on 20 year timescales (as in 
the dynamo, Kumar, Talon & Zahn 1999) or 1 year timescales (as in the 1.3 year oscillation, 
Kim & MacGregor 2001). We have reviewed this mechanism here so comparisons between 
this process and the zonal flow oscillations seen in the radiative region of our model can be 
clearly made. 



3. Numerical Model 

The numerical technique and model setup are identical to those in Rogers & Glatzmaier 
2005, except here we impose the equatorial rotation profile as a function of radius as inferred 
from helioseismology (fi( r ) * s specified to be 465nHz in the convection zone, 435nHz in the 
stable region and the tachocline is fit to an error function). We solve the Navier-Stokes 
equations with rotation in the anelastic approximation in 2D cylindrical geometry (r,(/>). 
The curl of the momentum equation, i.e., vorticity equation, is: 

duo _ . ^. . n dQ(r) g dT 1 dT dp „ 2 

The heat equation is: 



dT dT 

— + (v- V)T = -Ur(^r - (7 - l)Th p ) + (7 - l)Th p v r + 



dT — BT (1 

7 «[ V 2 T + (h p + h K ) — ] + 7 «[ V 2 T + (h p + h K ) — ] + f (4) 

In these equations, v is the velocity, with radial, v r , and longitudinal, vg, components. 
The vorticity is uj = V x v and is normal to the equatorial plane in this 2D geometry. The 
functions h p = dlnp/dr, h K = dlnR/dr, g (gravity), T (temperature), p (density) and 7 
(ratio of specific heats, c p /c v ) are radially dependent and taken from the solar model. T is 
the temperature perturbation and p is the pressure perturbation, which like u, are functions 
of r, 6 and time (t). 
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In our model we specify the thermal diffusivity as that given by the solar model, mul- 
tiplied by a constant for numerical stability: 

, 16(7 ^ 3 

k = kapmult * _ 2 _ — (5) 
3p JvCp 

here kapmult is generally set to 10 5 , a is the Stefan-Boltzman constant and k is the opacity. 
The viscous diffusivity, V = p/p, is set so that the dynamic viscosity (//) is constant. 

We calculate the pressure term in (3) using the longitudinal component of the momen- 
tum equation: 

1 dp dv e - -rm2^ h P dv r, 

These equations are supplemented by the continuity equation in the anelastic approximation 

V • pv = (7) 

which is satisfied by expressing pv as the curl of a streamfunction. The model extends from 
.001i? o to O.93i? . The subadiabaticity in the radiative region is given by the solar model 
and we specify the superadiabaticity in the convection zone to be 1CT 7 . 

These equations are solved using a Fourier spectral transform method in the longitudinal 
(0) direction and a finite difference scheme on a non-uniform grid in the radial (r) direction. 
Time advancing is done using the explicit Adams-Bashforth method for the nonlinear terms 
and an implicit Crank-Nicolson scheme for the linear terms. The boundaries are impermeable 
and stress-free. The inner boundary is isothermal and the outer boundary is held at a 
constant heat flux. This code is parallelized using message passing interface (MPI) and the 
resolution is 2048 longitudinal zones x 1500 radial zones, with 620 radial zones dedicated to 
the radiative region. In the region just below the convection zone the radial resolution is 
170km. This model was evolved for 1 simulated year, requiring nearly six million 5 second 
timesteps. 



4. Convection zone and Tachocline 

4.1. Angular Velocity Variations 

Helioseismic observations indicate that, in the equatorial plane, the convection zone 
spins faster than the radiative interior. The maintenance of differential rotation within the 
equatorial plane of the convection zone is a 3D process involving the transport of angular 
momentum in both radius and latitude and the Coriolis forces resulting from axisymmetric 



-7- 



meridional circulation. Since our 2D geometry captures only the transport in radius, we can 
not expect to achieve a realistic differential rotation profile with the 2D model. Therefore, 
we impose the observed equatorial angular velocity as a function of radius in our model. 
As mentioned above, the background angular velocity is set to a constant 465nHz in the 
convection zone and to a constant 435nHz in the radiation zone; we prescribe a smooth fit 
between these two values through the thin tachocline shear layer. The resulting gravity wave 
spectrum and angular momentum transport is then investigated. 

The time series of angular velocity, fi'(r, t), relative to the prescribed background an- 
gular velocity, £l(r), for this model is shown in Figure 1 over one simulated year. As seen 
in this figure, angular velocity in the convection zone is initially prograde, relative to the 
prescribed fl(r), in the lower part of the convection zone and retrograde motion in the upper 
part. However, after the initial period, angular velocity becomes very time dependent and 
is dominantly retrograde in the lower part of the convection zone and prograde in the upper 
part. This profile is expected, given the density stratification and rotation of the solar inte- 
rior (Glatzmaier et al. 2005). Intermittently prograde motion from the top of the convection 
zone will extend to the base of the convection zone and overshoot into the tachocline. The 
angular velocity variations produced by these motions vary in amplitude between +/- 15nHz, 
slightly larger than the amplitudes of the observed 1.3 year oscillation (Howe et al. 2001). 
Figure la clearly shows that the angular velocity of the tachocline mimics the behavior in 
the lower part of the convection zone; when the lower part is prograde, the tachocline is 
prograde and vice versa. This indicates that the behavior of the tachocline is dictated by 
convection zone dynamics, rather than by gravity waves (see below). 

In the overshoot region fluid motions are strongly nonlinear (Rogers & Glatzmaier 
2005b). Figure 2 shows the ratio of the horizontal fluid velocity to the horizontal phase 
speed (the Froude number) for a typical frequency and wavenumber (20/iHz, 1=10); this 
provides a measure of the nonlinearity of waves. For linear waves u x /c x « 1. However, 
as is clearly seen in the figure this criterion does not hold just below the convection zone. 
While only one ratio is shown as a function of radius, this can vary by an order of magnitude 
depending on the choice of frequency and wavenumber. 1 However, using reasonable values 
for frequency and wavenumber, the smallest value of this ratio just below the convection 
zone is 0.1, still far from meeting the linearity criterion above. Furthermore, it is clear in 
this figure that linearization may not be justified in the solar core either, as suggested by 
Press (1981). 



1 For larger (smaller) values of horizontal wave mode number this ratio is larger (smaller). Similarly, for 
larger (smaller) frequencies this ratio is smaller (larger). 
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4.2. Reynolds and Viscous Stresses 

The mean zonal flow is determined by a balance between viscous and Reynolds stresses 
(as seen in equation (1)). In order to better understand the angular velocity profile in 
time (as shown in Figure 1) we examine the sources of that balance. Figure 3 shows the 
horizontally averaged convergence of Reynolds stress (first term on the right hand side of 
(1), but now in cylindrical coordinates), the viscous stress (second term on the right hand 
side of (1), in cylindrical coordinates) and the sum of these (the time rate of change of the 
mean zonal flow U) in both the convection zone and overshoot region over five days. 

In the model's convection zone, the Reynolds stresses are much larger than the viscous 
stresses. The pattern of Reynolds stress convergence has a semi-periodic behavior with 
disturbances propagating up and down from a radius near l/?id cz [d cz is the convection 
zone depth). When positive (i.e., prograde) Reynolds stress convergence moves upward 
from this radius a negative value moves downward and vice versa. This process reverses 
with a frequency similar to the convective turnover frequency, roughly 10 times the inertial 
wave frequency. This oscillatory behavior is observed in our simulations only when we 
include a radiative region beneath the convection zone; we do not see it when we impose an 
impermeable lower boundary on the convection zone. Note also that this radius of l/3d cz is 
approximately where vg for the large eddies changes sign and below which smaller, counter- 
rotating cells exist. Figure 4 shows a zoomed-in region of the convection zone, displaying 
the vorticity within a few convective cells. Beneath the large cells spanning most of the 
convection zone is counter-rotating fluid that is driven by the plumes penetrating slightly 
into the stable overshoot region. As a plume rebounds it is diverted laterally and counter- 
rotating eddies are generated. It remains to be seen how these counter-rotating cells affect 
tachocline dynamics in a 3D simulation. 

Now consider the overshoot region in our model. Regions with positive (negative) 
Reynolds stress convergence in the lower convection zone correspond to the same signed 
disturbances in the overshoot region. These disturbances have a slight tendency to prop- 
agate downward, i.e., away from their source in the convection zone. When viscous terms 
are added (which in the overshoot region have magnitudes similar to the Reynolds stress 
terms) disturbances clearly propagate away from the convection zone (Fig. 3). This motion 
away from the convection zone 2 is also seen in Figure 1. The semi-periodic oscillations 
observed in dU/dt affects only the amplitude of U, but not the direction. That is, the 
semi-periodic oscillation seen dU / dt is not realized in U because of a slight asymmetry in 
the driving amplitude of prograde and retrograde motions. If, as has been assumed in the 



2 However, even if diffusion were neglected the disturbances do not move upward as in the QBO. 
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past (Kumar, Talon & Zahn 1999, Kim & MacGregor), these motions were produced with 
equivalent amplitudes, the semi-periodic oscillation observed in dU / dt would be reproduced 
in U, however, that is not the case here. 

Given the above results it is clear that overshooting plumes themselves transport angular 
momentum in a way which is distinct from the nonlinear interaction of low amplitude waves. 
The main properties of a wave driven oscillation (motion in time toward the source, timescales 
different than the forcing timescales, and a double peaked layer) are not seen here. Therefore, 
given the link between Reynolds stresses in the convection zone and those in the tachocline, 
it is clear that the overshooting plumes provide the main transfer of angular momentum to 
the tachocline. These plumes are buoyantly braked, radiatively damped and the nonlinearity 
allows for mode-mode transfer which affects the damping rate. Furthermore, the momentum 
transport by these plumes is time dependent and complex. Their effects need to be included 
when considering the dynamics of the tachocline and their influence should not be treated 
simply as an increased flux in linear waves. In fact, the linear approximation of waves is not 
justified in the overshoot region and tachocline. 



5. The Deep Interior and Core 

A banded radial differential rotation profile is observed deep within the radiative interior 
(Figure 5); however, the amplitude of the angular velocity variation is extremely low (0.02 
nHz) relative to what it is in the convection zone. This radial shear grows in time (Figure 6) 
due to wave-wave interactions which enforce shear. However, there is some indication that 
the growth is slowing and possibly reversing (Figure 6) due to increased viscous diffusion 
as angular velocity gradients increase. The banded differential rotation is due to the broad 
spectrum of frequencies generated and their continual deposition of angular momentum due 
to radiative diffusion as they propagate. Banded differential rotation in the sun's radiative 
interior has not been inferred from helioseismology, possibly because the amplitude is too 
low, or because these small amplitude disturbances are easily neutralized in 3D. In the core 
of our model the amplitude of the angular velocity variation is 2-3 orders of magnitude 
higher than in the bulk of the radiative interior. This indicates that angular momentum 
transport from waves to zonal flow is more efficient in the core, which could be due to 
a number of factors, including: wave breaking, critical layers (where the phase speed of 
the wave approaches the mean rotation rate) or inefficient reflection as the wave frequency 
approaches the Brunt- Vaisala frequency. Discerning which of these processes is dominant in 
this region is extremely difficult when a large spectrum of waves is continually being excited 
and is beyond the scope of this paper. 
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Maxima in the convergence of Reynolds stress deep within the radiative interior (but 
a significant distance below the overshoot region) propagate upward (Figure 7) as expected 
from the nonlinear interaction of low amplitude waves. In addition, significant interference 
is seen in the lower radiation zone, which likely comes from the interaction of inward prop- 
agating and outward reflected waves. While we have not observed a complete reversal of 
the angular velocity in the bulk of the radiative region, there is some evidence that one is 
underway (Figure 6) 3 . In the core the mean angular velocity switches from retrograde to 
prograde at nearly the same time as the tachocline switches from prograde to retrograde. 
This suggests that the angular velocity of the core is linked to the angular velocity of the 
tachocline and hence, of the convection zone. The link between tachocline angular velocity 
and that of the core due to selective filtering has been elucidated previously in Talon, Kumar 
& Zahn 2002. Here we make the link between the tachocline and the convection zone and 
propose that a self-consistent study of differential rotation in the sun must not treat the 
radiative and convective regions separately. 



6. Discussion 



Despite some similarity between the convective-radiative interface in the sun and the 
Earth's tropopause, there are several obvious differences. In the sun, convection is constantly 
driving gravity waves everywhere below the overshoot region; whereas the generation of large- 
amplitude gravity waves in the Earth's atmosphere is intermittent in time and space. Solar 
gravity waves travel down into a converging region of increasing density (and therefore, have 
very little chance of breaking); whereas in the Earth waves travel up into an expanding region 
of decreasing density. Furthermore, previous numerical simulations of convective penetration 
in a stratified atmosphere demonstrate a remarkable difference between penetration into an 
overlying stable region and penetration into an underlying stable region (Hurlburt, Toomre 
& Massaguer 1986). Penetration into an underlying stable region is characterized by thin 
localized downflows; whereas penetration into an overlying stable region is characterized by 
larger scale broad upflows. This asymmetry allows descending plumes to travel farther into 
an underlying stable region than ascending motions travel into an overlying stable region. 
Therefore, it is likely that penetrative convection plays a more crucial role at the base of 
the solar convection zone than it does at the Earth's tropopause. These differences can have 
profound effects on the role of overshoot and, hence, on the scale, frequency and amplitude 
of the waves generated and, so, on the angular momentum transport by these waves. 



3 It appears that in the deep interior a QBO-like oscillation is possible 
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Numerical simulations (Wedi & Smolarkiewicz 2005) of the Plumb-McEwan laboratory 
experiment attempting to reproduce the QBO, have shown that the type and period of an 
oscillation in the differential rotation profile depend sensitively on the forcing. In particular, 
it is found that random forcing rarely produces a periodic oscillation. Given the turbulent 
nature of the sun, it is likely that the forcing is fairly random. 

For reasons stated above, it is unlikely that there is a QBO-like oscillation associated 
with the solar tachcocline. However, QBO-like oscillations may occur in stars with radiative 
envelopes because of the inefficiency of overshoot into an overlying stable region, and because 
of more similar geometry. 

7. Conclusions 

We have presented self-consistent numerical simulations of convective overshoot and 
gravity wave generation and the angular momentum transport by these processes in a 2D 
model of the dynamics in the solar equatorial plane. We find that angular velocity variations 
in the tachocline are driven by angular momentum transported by overshooting plumes rather 
than nonlinear interaction of low amplitude waves. These overshooting plumes are strongly 
nonlinear disturbances, which can not be accurately represented as an increased flux of linear 
waves. We observe a semi-periodic oscillation in amplitude, but not in direction, of the mean 
flow in the tachocline because of an asymmetry in the driving of prograde and retrograde 
motions. Since we find that linear gravity waves are not dominant in the tachocline it is 
unlikely that they are responsible for the 1.3 year oscillation or the 11 year solar cycle. It is 
no surprise that overshooting motions play a dominant role in the tachocline and we expect 
these results will persist in three-dimensions. 

In the deep radiative interior the continual deposition of angular momentum by the 
nonlinear interaction of gravity waves produces a radially banded differential rotation. How- 
ever, it remains to be seen whether this pattern persists in 3D, considering it's very low 
amplitude. In the model's core, the amplitude of the differential rotation (i.e., angular ve- 
locity) is larger, about two orders of magnitude larger than that in the bulk of the radiative 
region and similar to the magnitude within the convection zone. We observe retrograde 
motion in the core reversing to prograde motion in step with the counter-reversal (prograde 
to retrograde) at the tachocline. When there is predominantly prograde flow at the base of 
the convection zone it selectively filters out prograde propagating gravity waves, allowing 
predominantly retrograde waves to propagate to the core where they deposit their (negative) 
angular momentum, and vice versa (Talon, Kumar & Zahn 2002). Like previous results our 
simulations therefore suggest that the angular velocity variations in the solar core are linked 
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to those in the tachocline. However, unlike previous results we show here that the variations 
in the tachocline are driven by convective overshooting, therefore, linking core rotation to 
convective motions. 
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Fig. 1. — Top: Angular velocity variations, relative to the prescribed solar profile, as a func- 
tion of time and radius. Red represents prograde motion, while blue represents retrograde 
motion. Black line represents the convective-radiative interface. Variations are shown in 
nHz. Bottom: Angular velocity variation as a function of time at the convective-radiative 
interface. Dotted line represents zero fluctuation about the mean. 
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Fig. 2. — The ratio of fluid velocities to the phase velocity (for a typical frequency and 
wavenumber), also known as the Froude number, a measure of the nonlinearity of the waves. 
This figure shows that motions are clearly nonlinear, both in the tachocline and the core. 
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Fig. 3. — Convergence of the Reynolds stress, Viscous stress and the sum, dU/dt, over a 
five-day period in both the convection zone (top) and overshoot region (bottom). A semi- 
periodic oscillation is seen both in the convection zone and overshoot region. This oscillation 
causes an oscillation in the amplitude (but not direction) of the mean zonal flow seen in Fig. 




Fig. 4. — Zoom in of vorticity in the convection zone. Large cells span the bulk of the con- 
vection zone with counter-rotating cells beneath due to the deflection of descending plumes 
as they encounter the stiff radiative region. 
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Fig. 5. — Angular Velocity variations in the deep interior. Banded radial differential rotation 
is seen. The deviations in the core are 2-3 orders of magnitude larger than those in the bulk 
of the radiative interior. There appears to be a link between deviations in the core and those 
in the tachocline (see Figure 1). 
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Fig. 6. — Time evolution of angular velocity variations. Some indication of reversal is seen 
in the bulk of the radiative interior between 223 days and 348 days. In the core there was 
clearly a reversal, which Figure 5 tells us was around 130 days. 
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Fig. 7. — Convergence of the Reynolds stress in the bulk of the radiative interior (blue 
represents negative, red positive values). Disturbances move upward in time as expected 
from the nonlinear interaction of low amplitude waves. Significant interference is seen at 
smaller radii, probably due to the interaction of inward propagating and reflected waves. 



